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PSEUDO WIGNER-VILLE DISTRIBUTION, 
COMPUTER PROGRAM AND ITS APPLICATIONS 
TO TIME-FREQUENCY DOMAIN PROBLEMS 

by 

Jae-Jin Jeon and Young S. Shin 



ABSTRCT 

Machinery operating in non-stationary mode generates a 
signature which at each instant of time has a distinct frequency. 
A time-frequency domain representation is needed to 
characterize such signature. Pseudo Wigner-Ville distribution is 
ideally suited for portraying non-stationary signal in the time- 
frequency domain and carried out by adapting the fast Fourier 
transform algorithm. The important parameters affecting the 
pseudo Wigner-Ville distribution are discussed and sensitivity 
analyses are also performed. Practical examples of an actual 
transient signal are used to illustrate its dynamic features 
jointly in time and frequency. 



I. INTRODUCTION 

The physical condition or state of health of machineries which 
operate in transient or other non-stationary modes are difficult to predict 
with any degree of accuracy. It is common to practice periodic preventive 
maintenance on these machineries in order to avoid failures and prolong the 
useful operating life of the equipment. 
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In order to assess the physical condition of machinery without 
complete disassembly, a physical measurement of its vibrations is conducted 
using an accelerometer. Other sensors, such as temperature or pressure 
transducers, could also be used. There are other methods, including motor 
current signature analysis on electrically driven machinery and wear debris 
analysis which could be used. However, vibrations are used predominantly 
for machinery condition monitoring. The vibrations are recorded in the time 
domain. 



There is a need for a method to represent the time dependent events 
which occur with machinery operating in non-stationary modes. At each 
instant in time as the speed of the machinery changes, the frequency content 
will also chamge. The pseudo Wigner-Ville distribution(PWVD) is the method 
which was chosen to portray these time dependent changes. This is a 
continuation of work initially performed and published by Rossano, Hamilton 
and Shin [1]. 

The pseudo Wigner-Ville distribution is a three dimensional (time, 
frequency, amplitude) representation of an input signal and is ideally suited 
for describing transient or other non-stationary phenomena. The Wigner 
Distribution (WDF) has been used in the areas of optics [2,3,4] and speech 
analysis [5,6]. Wahl and Bolton [7] used it to identify structure-bome noise 
components. Flandrin et. al. [8] recently proposed its use in the area of 
machinery condition monitoring and diagnostics, while Forrester [9] is 
investigating its use in gear fault detection. 

For such a non-stationary signal analysis, spectrogram is commonly 
used, which is based on the assumption that it is a collection of a short 
duration stationary signals. A major drawback of this approach is that the 
frequency resolution is directly affected by the duration of short stationary 
time, which subsequently determines the time resolution. A method for time- 
frequency domain signal characterization that overcomes this drawback is 
the Wigner distribution which was first introduced by Wigner [10] in 1932 to 
study the problem of statistical equilibrium in quantum mechanics. The 
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frequency and time resolutions of the Wigner distribution are not determined 
by the short duration but rather determined by the selection of desired 
resolution of the signal itself. 

This paper discusses the important parameters affecting the PWVD 
in order to machinery condition monitoring and presents numerical examples 
of PWVD using synthetically generated signals. It is found that the PWVD 
is very effective in machinery condition monitoring which operates in non- 
stationary modes. 
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II. PSEUDO WIGNER-VILLE DISTRIBUTION 



A. Wigner Distribution Function 

Signal associated with most vibrational phenomena are in general time 
varying, which means that their characteristics change with time and they 
have various features in different time frames. The general spectrogram 
usually requires a large time-bandwidth product to reduce the estimated bias 
and variability. In the case of a signal containing some transients or 
nonstationary conditions, the traditional approach in signal analysis fails to 
describe the dynamics of the signal's frequency component changes. 

The general expression of the time-frequency distribution of a signal, 
w(t,(0) is given by, [11] 

w(t,co) = — s*(u -— ) s(u-i-— ) du dx d0 (1) 

2n 2 2 

where s(u) is the time signal, s*(u) is its complex conjugate, and (f>(6, t' is an 

arbitrary function called the kernel. By choosing different kernels, different 
distributions are obtained. Wigner distribution is obtained by taking (j)(9,T) 

= 1. The range of all integrations is from - oo to «« unless otherwise noted. 
Substituting the kernel <^(0, t) = 1 to Eq. (1), 

w(t,(o) = jjs*(u-^) e~-^^^5(u-t) s(u-i--^) dxdu (2) 

From Eq. (2) the Wigner distribution is obtained, 

w(t,CD) = dx (3) 
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One of the basic frequency representations of a signal is the power 
density spectrum, which characterizes the signal’s frequency component 
distribution. The power spectral density function S(ty) of a signal s(t) can be 
related to the Fourier transform of the signal's autocorrelation function R(t): 



S(co)= dr 

with 

/?(t) = Js(t) s(t + T)dt . 



(4) 

(5) 



From this relation a time-dependent power spectral density function 
can be written as 



w(t,( 0 ) = jRt(x) e j“^dx (6) 

where now R^(z) is a time-dependent or local autocorrelation function. Mark 
[12] argued for symmetry, 

fi,(T) = s'(t-|)s(t + |) (7) 

which gives the Wigner distribution function. 



B. Properties of the Wigner distribution function 

The properties of the WDF [13,14] are summarized and reinterpreted 
with this new formulation as follows; 

(i) The WDF is a real-valued function. 

w*(t,tw) = w(t,0)) (8) 
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(ii) The integral of the WDF with respect to frequency and time yields 
the instantaneous signal power and the signal's power density spectrum, 
respectively. 

o= 2 

|w(t,co)dt» = 2;r|s(t)| , 

— OO 

2 (9) 

oo ^ 

Jw(t,fti)dt = 2;r|S(0(i)| . 

-oo 



(iii) A time or frequency shift in the signal have the shift in the WDF. 



If s(t) — > s(t + to), then w(t,to) w(t + to,ty), (10) 

if s(t) then w(t,ft)) w(t,0)+(Mo)- (ID 

(iv) The WDF is symmetrical in time for a given signal. 

If s(t) — > s(-t), then w(t,to) w(-t,Ct)), (12) 

if s(t) — > s *(t), then w(t,ft)) — > w(t,-ft)). (13) 



(v) The WDF is not always positive. 



(vi) The integration of the square of the WDF equals the square of the 
time integration of the signal's power. This is the counterpart of Parseval's 
relation of the WDF, called Mayol's fomula. 



oo ^ 

J|w(t,<»)| dtdCD = 

— oo 



Js’COdt 



(14) 



(vii) The WDF possesses basic symmetry under the interchange of time 
8ind frequency parameters with the Fourier transform of a given signal. 

Let 



s(t) = — Te-’"' S(6i) dto; 
In ^ 



(15) 
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then 



w(t,to) 



— ? S(to + -) S*((y - -) d^. 
iK^ 2 2 ^ 



( 16 ) 
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III. IMPLEMENTATION WITH DIGITAL SIGNAL PROCESSING 



A. Discrete Wigner Distribution Function 

There are two distinct advantages for the calculation of the WDF. 
First, it has the form of the Fourier transform and the existing FFT 
algorithm can be adapted for its computation. Second, for a finite time signal, 
its integration is finite within the record length of the existing signal. 

The discrete time Wigner distribution as developed by Claasen and 
Mecklenbrauker [13] is expressed by, 

w(t,co) = 2 ^ire~j^“Vt + x)s*(t-x) (17) 

X=-oo 

The discrete version of Eq. (17) for a sampled signal s(n), where n=0 to N-1, 
has the form, 

1 N-1 ^ -j— nk 

w(f,k) = — £ s(^ + n) s (^-n) e N ^ k=0,l,2,...N-l (18) 

N„=o 

where s(m)=0 for m < 0 and m > N-1. However, in order to utilize the FFT 
algorithm, it must be assumed that the local autocorrelation function has a 
periodicity of N. This is just for operational convenience and should not apply 
to the interpretation of s(m). Eq. (18) can be rewritten as, 

1 N'l * -y— n(k+m— ) 

w[^,k -I- m(N / 2)] = — Ts(.^-i-n)s (f-n)e ^ ^ 

Nn=o 

1 N-1 , -j— nk . . 

= -Xs(^ + n)s (^-n)e N (19) 

N n=0 

= W(f,k) 
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since = 1 for m=integers. 

Eq. (19) indicates that the WDF has a periodicity of N/2. Hence, even 
when the sampling of s(t) satisfies the Nyquist criteria, there are still aliasing 
components in the WDF. A simple approach to avoid aliasing is to use an 
analytic signal before computing the WDF. In 1948, J. Ville [15] proposed the 
use of the analytic signal in time-frequency representations of a real signal. 



B. Analytic Signal 

An analytic signal is a complex signal which contains both real and 
imaginary components. The advantage of using the analytic signal is that in 
the frequency domain the amplitude of negative frequency components are 
zero. This satisfies mathematical completeness of the problem by accounting 
for all frequencies, yet does not limit the practical application since only 
positive frequency components have a practical interpretation. The 
imaginary part is obtained by Hilbert transform. The analytic signal may be 
expressed by, 

S(t) = Sr(t) -I- j H{Sr(t)} (20) 

where H{Sj-(t)} is a Hilbert transform and generated by the convolution of 
the impulse response h(t) of a 90-degree phase shift as follows: 

H{s,(t)} = s,(t) * h(t) 

h(t) = ^ sin^(m/2) 

;rt 

= 0 , 

where * denotes the convolution. Rewriting Eq. (21) to discrete form, 

oo 

H{Sr(n)} = Xh(n-m) Sr(m) (22) 

m=-oo 



( 21 ) 

t = 0, 
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The distribution resulting from an analytic signal being processed through 
the Wigner distribution is commonly termed as Wigner-Ville distribution. 



C. WDF with Digtal Signal Processing 

To calculate the Wigner distribution of the sampled data, it is 
necessary that Eq. (18) be modified to Eq. (23), because the WDF has N/2 
periodicity. 

2N-1 

w(mAt,kAo)) = 2At ^s[(m + n)At] s*[(m-n)At] (23) 

n=0 

where Aco = 7C / (2NAt) and At is the sampling interval. The algorithm used 
in this paper is based on one written by Wahl and Bolton[7] and can be 
expressed as: 

w(mAt,kA6)) = Re [2At FFT(corr(i))] 

coir(i) = s(m + i-l) s*(m-i + l), m > i ( 24 ) 

= 0, m < i 

where 1 < i < N + 1 

coit(2N - i + 2) = corr*(i), 2<i<N 

The frequency resolution, Aco , in Eq. (23) is different from that obtained by 
FFT of the original N point time record in two respects. The first difference is 
that the argument of the time signal and its conjugate contains a factor of 
1/2, and secondly, the autocorrelation of the time signal is twice the length of 
the original record and therfore the FFT is evaluated over 2N points. The 
result is, that the WDF frequency resolution is one forth the resolution of an 
ordinary power spectrum density fimction. 

Before processing the WDF, a modified Hamming window is applied to 
the time domain signal to reduce the leakage caused by the discontinuity of 
the finite record of data, which will be called as data tapering. This type of 
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window is preferable since it alters the amplitude of fewer data points at the 
beginning and the end of the data block. A modified Hamming window, D(t) is 
given by: 



0.54 - 0.46 * cos(lOrtbT), 
D(t) = { 1.0, 

0.54 - 0.46 * cosdOrt(T-t)yT), 



0 < t < T/10, 

T/10 < t < 9T/10, (25) 

9T/10 < t < T. 



D(t) 




t 



Fig. 1 Modified Hamming Window 



Two other characteristics of the WDF should be also noted. First, the 
WDF of the sum of two signals is equal to the sum of the WDF of each signal 
plus cross term that appear when the cross-correlation of the two signal is 
non-zero. Second, the WDF may have negative values, which may be largely 
caused by interference due to the presence of these cross terms. In the case 
of input signals that contain multi-frequency components, the Wigner-Ville 
distribution of most signals are very complicated and difficult to interpret. 

There are two methods to suppress the interference components of the 
WDF. Claasen and Mecklenbrauker[12] describe the application of a sliding 
window in the time domain before calculating WDF. The WDF obtained vvdth 
a window function is called the Pseudo- Wigner distribution function. A 
second option is to smooth the WDF with a sliding averaging window in time- 
frequency plane. In both case the result is to deemphasize components 
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arising from calculations and to emphasize deterministic components. 
Obviously, averaging a Wigner-Ville distribution will result in a Pseudo 
Wigner-Ville distribution. 



In this research, a sliding exponential window in the time-frequency 
domain was chosen. That is, a Gaussian window function, G(t, co) is selected 
to reduce the interference and to avoid the negative values as follows: 
let 



G(t,(0) 



then 



1 



-( 



0 ) 



2af 2a: 






0) 



(26) 



w(t,(o) = — (o') G(t-t',(o-(o') dt'dco' 

2k 



>0 



(27) 



where O^, > 0 and 1/2 [16]. The time and the frequency 

resolution At and A(0 of this Gaussian window are related by, 

Oj=jAt, Oqj = k Ao) (28) 

in the discrete form. Then the condition for the WDF to be positive in this 
case is 

j At k Aco > 1/2. (29) 

This is the time-frequency version of Heisenberg's uncertainty relation[14]. 
If the segmentation of time and frequency for a given signal from Eq. (3) 
violates this uncertainty principle, the corresponding WDF may not be 
positive. 

To perform the convolution on the sampled WDF, the Gaussian window 
function was applied to the range ±20{ and ±2(7 qj. Selecting co and t to be 

the multiple of time and frequency steps, the sampled Gaussian window 
function is expressed by. 
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(30) 



G(p,q) = 



1 

2 tc j k AtAco 



rz 

•2 



.2 A 



2j 



2k^ 



where p and q are an integer numbers in the range ±2j and ±2k, respectively. 
The convolution of the sampled WDF and the Gaussian window function can 
be evaluated as follows; 

At Aco 

w'(tm) = — - — X Xw(p,q) G(p-^, q-m) (31) 

271 p=f-j q=m-k 

where w'(^,m) is the smoothed WDF or Pseudo Wigner-Ville distribution. 



Fig. 2 shows a block diagram for computational sequence of the Pseudo 
Wigner-Ville distribution. A time-varying signal sampled with the Nyquist 
rate is first high passed through a digital filter if the signal involves the zero 
frequency component, i.e., DC component, and converted into the analytic 
signal through a Hilbert transform. Then, the time-dependent correlation 
function is computed and the result is the WDF in terms of both time and 
frequency domain by FFT. The final step is to compute the convolution with a 
Gaussian window. 




Fig. 2 Computational block diagram of Pseudo Wigner-Ville Distribution 
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D. Highpass Digital Filter 



Filters are a paticularly important class of linear time-invariant 
system. Strictly speaking, the term frequency-selective filter suggests a 
system that passes certain frequency components and totally rejects all 
others, but in a broader context any system that modifies cetain frequencies 
relative to others is also call filter. The design of filters involves the following 
stages ; (1) the specification of the desired properties of the system; (2) the 
approximation of the specifications using a causal discrete-time system; and 
(3) the realization of the system. 



In this paper, Nonrecursive(finite impulse response - FIR) highpass 
filter was used for the elimination of undesired low frequency components. 
The basic design was to use a symmetric filter of the form. 



M 





y(i) = X s(i-k) 


(32) 




k=-M 




with 


t>-k = 


(33) 


and 


, sin 2;rBkT 

bk = 

7 ± 


(34) 



where is the filter weights, y(i) is the filtered signal, B is the cutoff 
frequency, T is a sampling interval and M is the span of the filter; 2M-I-1 
weights are employed because of symmetry, only M-t-1 need be generated. 

The weights are computed over the range -M to M. The weights are 

multiplied by a window function. Potter discusses a number of windows in 
the referenced work. His P310 window was found to be appropriate for filter 
implementation. It takes the form. 






where 




3 

do + 2 X dp cos 

p=-3 



;rpk 



(35) 
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and 



c, . = - k = ± M 

^ 2 

= 1 otherwise 

do = 1 

d. i = di = 0.684988 

d.2 = d2 = 0.202701 

d.3 = dj = 0.0177127 



3 

w = do + 2Xd = 2.8108034 

p=-3 



( 36 ) 



(37) 



For a highpass filter with pass band from the cutoff frequency(B) to the 
maximum frequency, generate a low pass filter on the range 0 - B, and then 
subtract the central weight from unity and change the signs of the remainder 
of the weights. 
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IV. EXAMPLES AND DISCUSSIONS 



Machinery operating in transient mode generates a signature in which 
the frequency content varies at each instant of time. To characterize such 
signatures and to understand the vibrational behavior of such machineries, 
time-frequency domain representation of the signal is needed. As discussed 
in the previous sections, Wigner distribution is a signal transformation that 
is particularly suited for the time-frequency analysis of nonstationary signals. 
There are many advantages of using PWVD for both steady and transient 
signals. However, there are also several disadvantages, for example, the 
drastic increase of peak value when the frequency content of signal changes 
abruptly. A computer program has been developed for PWVD. Two different 
versions are available at the present time; workstation and IBM PC 
compatible. 



A. Harmonic Wave 

Fig. 3 shows the PWVD of the pure sine wave with two frequency 
components (lOOHz, 400Hz), respectively. The modified Hamming window 
was applied to the time domain signal and the Gaussian smoothing window 
function was applied on time-frequency domain Winger-Ville distribution. 
The slope of the end edges are due to data tapering by using the modified 
Hamming window. Fig. 4 shows the PWVD of the sine wave that have the 10 
% and 50 % signal to noise ratio, respectively. The shape of PWVD is 
changed at the crest by the contamination of noise. The crest has the 
complicated shape with decreasing of signal to noise ratio. However, the ' 
PWVD well represents the signal components from the given signal with 
noise. The practical example is shown in Fig. 14. The notation fg and N 
used in the figures are sampling frequency and the total number of time data 
points. 
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«Q 




Fig. 3 Pseudo Wigner-Ville distribution of 100 and 400 Hz Pure Sine Waves 
(f =1000 Hz, N=256 and Smoothing Window Size=10xl0) 
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0.0100 



(a) 





Fig. 4 Pseudo Wigner-Ville distribution of 300 and 750 Hz sine waves; signal 
to noise ratio (a) 10 % and (b) 50 % . 

(fs = 2048 Hz, N=1024 and smoothing window size=18xl8) 
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B. Harmonic Wave with Stepwise Frequency Changes 



Fig. 5 shows the PWVD of the sine wave with 500 Hz in the time from 
0.085 sec to 0.17 sec. The PWVD well represents the time delay of the signal. 
The Fig. 6 shows (a) the sine wave with stepwise frequency changes, 100 Hz, 
250 Hz and 500 Hz and (b) its PWVD. The PWVD shows the time delay and 
frequency component of the signal. The wide spread of PWVD at the edge of 
each frequency region is noticed. This phenomenon is caused by the 
discontinuity of the signal in time domain and the leakage in digital signal 
processing. This effect may be reduced by applying the data tapering to the 
actual signal block. Nevertheless the PWVD represented the characteristics 
of the signal well. PWVD can portray the characteristics of the steady state 
signals involving time delay and multi-frequency components . If different 
size of the smoothing window are applied, the PWVD amplitude changes, but 
the total energy remains unchanged. 




Fig. 5 500 Hz sine wave with finite duration 
(fg=2000 Hz, N=512 and Smoothing Window Size=10xl0) 
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Magnitude 



1.5 
1.0 
0.5 
0.0 
- 0.5 
- 1.0 
- 1.5 

0.0 32.0 64.0 96.0 126.0 160.0 192.0 224.0 256.0 

Time (msec) 

(a) Time signal 



<0 




(b)PWVD 




T 



Fig. 6 Sine Wave with Stepwise Frequency Changes: 100, 250 and 500 Hz 
(fg=2000 Hz, N=512 and Smoothing Window Size=10xl0) 
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C. Composite Signal with Two Frequency Components at Each Time 

The PWVDs of the nonstationary signals were studied and the results 
were shown in Fig. 7 through 10. Fig. 7 shows (a) the time signal composed 
of two sweeping frequency components at each time, one increasing and the 
other decreasing with the same rate, and (b) its Wigner-Ville distribution 
(before applying the smoothing window) (c) its contour plot of WDF and (d) its 
pseudo Wigner-Ville distribution (after applying the smoothing window), 
respectively. 

The effect of cross (or interference) term is significant and appeared in 
the average frequency region. This is one of the disadvantages of using 
Wigner-Ville distribution but it is a characteristic of the distribution. When 
Gaussian window was applied to Wigner-Ville distribution, the effect of cross 
term disappeared. The main lobe of PWVD is wider and its amplitude is 
significantly reduced. The large peak at the intersection point of two 
sweeping frequency signals is mainly caused by the doubling effect of 
amplitudes of two signals. 




(a) Time signal 
Fig.7 (continued) 
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Utade 





(c) Contour plot of WDF 
Fig, 7 (continued) 
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(d)PWVD 



Fig. 7 Composite Signal with Two Frequency Components at Each Time 

s(t)=4cos(2;t 32t^) + 4 cos{2;t(40+32(2-t)]t} 

(f =256 Hz, N=256 and Smoothing Window Size=10xl0) 
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D. Linear Chirp Signal 



Another type of a non-stationary signal sweeps up and down in 
frequency is called a linear chirp signal and is shown in Fig. 8(a). This signal 
has only one frequency component at each time. The effect of cross terms 
appears in the Wigner-Ville distribution, as shown in Fig. 7(b). The 
smoothing window was applied to Wigner-Ville distribution and the result is 
shown in Fig. 7(d). Fig. 8 (c) is the contour plot of WDF. As expected, the 
effect of cross term is significantly reduced. However, the unusual peak 
(called 'ghost' peak) appeared at the point where the direction of sweep 
changes. To understand the cause of this phenomenon, the PWVD was 
integrated along the frequency axis and it was found that the square root of 
the resultant amplitude was the amplitude of original time signal, implying 
that the energy content remained constant. The following function was used 
to generate the linear chirp signal: 



s(t)= sin 2;rj^30 + , 

220(512 -i) 



s(t) = - sin 



2;d 30 + 



256 



(0.256 -t) 



1 < i < 256 
256 < i < 512 



(38) 



where t = (i-1) dt and dt=0.0005. 




Fig.8 (continued) 
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(b)WDF 




(c) Contour plot of WDF 
Fig.8 (continued) 
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(d)PWVD 



Fig. 8 Linear Chirp Signal with One Frequency Component at Each Time 
(f =2000 Hz, N=512 and Smoothing Window Size=16xl6) 



26 



E. Composite Signal of Sweeping-up and Steady Frequency 

The signal which sweeps up along the frequency for first 0.5 second 
and holds to a constant frequency for next 0.5 second was considered. This 
signal is typical speed profile of start-up stage of pump. Fig, 9 shows (a) 
PWVD and (b) its contour plot. The interesting phenomenon was observed in 
PWVD that the sweep-up portion of signal (first half seconds) has a lower 
amplitude and wider main lobe compared with the steady frequency region of 
signal (second half seconds). When the PWVD was integrated along the 
frequency axis and it was found that the resultant amplitudes in these two 
regions are same. The following functions were used to generate the desired 
signal: 



s(t) = 4cos(27c32t“), 0<t<0.5sec. 

s(t) = 4cos(27u64t), 0.5 < t < 1.0 sec. 

Fig. 10 is the PWVD of the signal which sweeps up along the fi*equency 
with a logarithmic rate, that is, the sweep rate is propotional to the square 
root of time. It was found that the maximum magnitude of the PWVD 
increases with increasing the frequency. This fact is shown that the PWVD 
of the stable signal has a larger magnitude than the unstable signals 
although having the same magnitude in time domain and the PWVD is the 
good tool for the analysis of the stability of the signal. The following 
functions were used to generate the desired signal: 

s(t) = 4cos{2;r(30 + 60 t*''^)t}. (40) 

The instantaneous frequency of the Eq.(40) is the derivative of the argument 
of the cosine function which is 

fit) = 90 4- 30 (41) 



27 





(b) Contour plot of PWVD 



Fig. 9 PWVD of a Composite Signal of Sweeping-up and Steady Frequency 
(f =256 Hz, N=256 and Smoothing Window Size=10xl0) 

o 
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0.4 




Fig. 10 PWVD of a signal of sweeping-up with a logarithmic rate with time. 
(fs=256 Hz, N=256 and Smoothig Window Size = 10x10) 
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F. Sweep Rate Efifect 



The effect of sweep rate on PWVD was investigated. The sweep rate is 
the frequency change per unit time. The power spectrum density of a typical 
swept sine is shown in Fig. 11. The incorrect assmnption is often made that a 
swept sine of constant amplitude has a flat spectrum. As can be seen from 
the plot this is not so. Fig. 12 shows the PWVDs of the linear chirp signal 
with a various sweep rates:(a) has zero sweep rate and (b) has lower sweep 
rate than (c). It can be seen that the amplitude of PWVD decreases with 
increasing sweep rate but energy remains unchanged. This result appeared to 
be caused by Heisenberg’s imcertainty relation between time and frequency. 
However, based on this study, it is clear that the 'ghost' peak (see Fig. 8) 
appears due to the instantaneous zero sweep rate at the point where the 
direction of sweep changes. Also the peak value is affected by the size of 
smoothing window. 




Fig. 11 Swept sine wave spectrum. 
(Frequency range 1 to 20 Hz, sweep time 22.5 sec) 
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mpliwde 




(c) s(t) = cos (2ji 64 t2) 



Fig. 12. The Effect of Sweep Rates To Pseudo Wigner-Ville Distribution 
(f =256 Hz, N=256 and Smoothing Window Size=10xl0) 
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G. Harmonic Wave with Some Glitches 



The interesting phenomena on the signal with an abnormal 
components as a fault were investigated. Fig 13 shows the PWVD of the 
harmonic wave with glitch at a small region of the time record: (a) is the time 
signal, (b) is the PWVD and (c) is the contour plot of the PWVD. It can be 
seen that the PWVD of the signal Fig. 13(a) well represents the loacation of 
each glitch and its frequency components. This characteristic of the PWVD is 
useful to detect the faults or glitch and to monitor the condition on the 
vibrational machinery having the periodicity such as a gear train. The 
general rotating machinery has a periodic signal pattern in time domain. 




Fig. 13 (continued) 
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Fig. 13 PWVD od the signal with gltches; (a) time signal, (b) PWVD and (c) 
contour plot of PWVD. 

(fs=256 Hz, N=256, smoothing window size 10x10) 



34 



H. Actual Fan Signal 



The acceleration signal of a fan was measured at the steady state 
speed and the result was shown in Fig. 14. The crest has the complicated 
shape on time axis as shown in Fig. 4. The first peak is the fundamental 
frequency of the blade rate and the second peak is 3rd harmonics. The third 
peak is the fundamental frequency of motor by the pole. The measured 
vibration signal was contaminated with the noise. If the measured signal 
involves the faults, the PWVD will represent the different pattern having the 
abnormal frequency components in comparison with the normal condition 
with time. 



'‘•"'00 „ 

1 




Fig. 14 PWVD of the actual fan signal in the constant speed, 
(fs = 512 Hz, N=512, smoothing window size 13x13) 
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I. Actual Pump Start-up RPM Signal 

The start-up transient speed of the pump was measured and the 
results were shown in Fig. 15. The time signal is measured by magnetic 
sensor and has the same magnitude independently on time. The PWVD is 
shown in Fig. 15(a) and the contour view is shown in Fig. 15(b). The contour 
plot shows that the speed of the pump runs up when initially started, reaches 
the maximum RPM and coasts down gradually. Near the maximum speed 
during the run up, the sweep rate was rapidly decreased and, as a result, the 
peak value was rapidly increased. When the sweep rate is close to zero at the 
normalized time of 0.4, the amplitude attains the maximum value. 



(a) 




Fig. 15 (continued) 
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Fig. 15. Pseudo Wigner-Ville Distribution of Transient Speed of the Pump; 
(a) PWVD and (b) contour plot. 



V. CONCLUSIONS 



The pseudo Wigner-Ville distribution has been investigated and applied 
to analyzing non-stationary signals typical of transient machinery signatures. 
The results of this research will be a valuable asset for condition monitoring 
of transient machinery. The following conclusions can be drawn: 

(1) The pseudo Wigner-Ville distribution is ideally suited for portraying non- 
stationary time signals. 

(2) The use of modified Hamming window to time signals is effective to 
reduce the edge effect of discontinuity. 

(3) The use of the analytic signal in calculating the Wigner distribution 
eliminates aliasing problem. 

(4) The Gaussian window function for smoothing the Wigner-Ville 
distribution is very effective and the presence of cross terms is significantly 
reduced. 

(5) Both the amplitude and the main lobe of the pseudo Wigner-Ville 
distribution is significantly affected by the sweep rate. As the absolute sweep 
rate increases, the amplitude of the PWVD decreases and the main lobe 
becomes wider. 

(6) The PWVD characterizes the time-frequency domain distribution of the 
signal well and may be useful tool for the machinery condition monitoring. 
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APPENDIX A. USER'S GUIDE OF THE PROGRAM 



This program calculates either the Wigner-Ville distribution or the 
smoothed Wigner-Ville(Pseudo Wigner-Ville) distribution(PWVD) of a time 
series. This program uses the FORTRAN 77 language and it is possible to 
run at workstation and IBM-PC compatible computer. The user supplies 
the real data of the time history, the sampling parameters, the digital filter 
parameters, the smoothing parameters and the output parameters. The 
program outputs 2-D array containing the Wiger- Ville distribution 
(rwdf.out) and the smoothed Wigner-Ville distribution, i.e, the pseudo 
Wigner-Ville distribution (rswdf.out) dependently output parameters. 



A-1. PROGRAM INPUT/OUTPUT (Main Program) 

This program is interactive and displays an explanations about a 
desired values. Fig. A-1 is a flowchart of this program and the 
computational block diagram is shown in Fig. 2. 

1) INPUT 

This program displays the input variables as follow. 

Enter name of signal input file 

The user must supply the input filename under 25 characters. 
Number of the sampled data point 

The user must input the number of sample data(dp). If dp is larger 
than 2048, the user must change np in a available memory size of the user's 
computer and recompile the source program. The dp must be a multiple of 
2 because of FFT. 



START ^ 




Fig. A-1 Flow chart (continued) 
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Fig. A-1 Flow chart 



The next step is, 



Do you wish to remove the mean value ? 

Enter 1 for yes or 0 for no. 

The user inputs the desired value in according to the characteristics 
of input time signal. If the time signal have the DC components, the user 
had better select mvopt=l. 



Do you want to apply highpass digital 
filter to the original data ? (y/n) 



If y' is selected, the desired cutoff (half-power point) frequency is inputed 
after displaying the follow command. 

Enter the cutoff frequency of 
the digital highpass filter (in Hz) 



This program uses the non-recursive(finite impulse response - FIR) 
highpass filter for the elimination of undesired low frequency components. 
If the user wants to pass the highpass filter, input the desired cutoff 
frequency in Hz. 

The next procedure is. 

Input the desired reduction size 
input 1 for 64 by 32 
input 2 for 128 by 64 
input 3 for 128 by 128 
input 4 for 256 by 128 
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This input statement describes the desired output file size for the output 
and the plotting because the original Wigner distribution involving the all 
domains needs the huge memory size. 

As the next step, this program chooses the 

The desired smoothing window size 

input the smoothing parameter for 
frequency in Gaussian function(nO 

input the smoothing parameter for 
time in Gaussian function(mt) 



The smoothing parameters for frequency (nf) and time (mt) depend on the 
number of the sampled data. The user inputs the integer values satisfying 
the following condition, 

nf * mt > — (A-1) 

n 

We recommend to select the square smoothing window size. 



2) OUTPUT 

This program generates the two output files. One is the file for the 
unsmoothed Wigner-Ville distribution and the other is the file for the 
smoothed Wigner-Ville distribution. The output files consist of one coliimn 
along with the time and frequency axis dependently the output parameter 
as following descriptions. 



O.OOOOOE-t-OO 

O.lOOOOE-t-Ol 



: initial time 
: time record length 



O.OOOOOE+00 : starting frequency 

0.25600E+03 ; end frequency 

128 64 : reduction size 

0.53830E-06 

0.49456E-06 

0.18271E-06 

: the results of PWVD or WVD 



The graphic output of the results used of CA-DISSPLA version 11.0. 
The 3-D graphic program describes the unsmoothed or smoothed Wigner- 
Ville distributionCor pseudo Wigner-Ville distribution). The 2-D graphic 
program describes the contour plot of the results and the time series. 
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A- 2. EXPLANATIONS OF INPUT VARIABLES 
Variable names Descriptions 



inname 


Input file containing the real time signal. 


dp 

mt 

nf 


The number of the sampled data point. 
Smoothing parameter for time in Gaussian Fn. 
Smoothing parameter for frequency in Gaussian 
Fn. 


mvopt 


Option with respect to removing mean value 
if mvopt = 1, then zero mean, 
if mvopt = 0, then no. 


bw 


Cutoff frequency of highpass filter 
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A-3. SUBROUTINES 



1) indata(dp,tin,ain) 

This subroutine reads the input file 'inname'. The input data file 
involves the time and amplitude in real value. The reading format is a free 
format. 



2) dtcalc(dp,tin,dt) 

This subroutine calculates the mean time interval dt from the input 
time signal. Because the obtained data from A/D converter or signal 
processor doesn’t have the exactly same interval. Delta time(dt) is given as 
follow, 

dt = (total record length) / (number of data point-1) (A-2) 



3) mean(dp,ain) 

This subroutine calculates the mean value and removes the mean 
value of the signal under the condition of the input variable mvopt. In the 
case the signal has the DC component, it is recommended to remove the 
mean value. A highpass digital filter must use for the elimination of DC 
component. 



4) filter(m,6,t,bk) 

This subroutine generates a nonrecursive (finite impulse response - 
FIR) filter weights and uses the method devised by Potter, Bickford and 
Glaze. Nonrecursive highpass filter was used for the elimination of 
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undesired low frequency components. The variables m, b, t, and bk are the 
number of the weights, the cutoff frequency of the filter, the time interval 
and the storage array of the filter weights, respectively. 



The basic design was to use a symmetric filter of the form. 



m 





y(i) = X bk s(i-k) 


(A-3) 




k=-m 




with 


b-k = bk 


(A-4) 


and 


, sin 2;r^kt 






bk = 


(A-5) 



where b^. is the filter weights, y(i) is the filtered signal, s(i) is the original 
signal, b is the cutoff frequency, t is a sampling interval and mo is the span 
of the filter; 2m+l weights are employed because of symmetry, only m+1 
need be generated. The weights are computed over the range -m to m. 

The weights are multiplied by a window function. Potter discusses a 
number of windows in the referenced work. His P310 window was found to 
be appropriate for filter implementation. It takes the form. 



w 



k 



where 




+ 



2 



3 

X dp cos 

p=-3 



;rpk 

m 




k = ± m 
otherwise 



do = 1 

d.i = di = 0.684988 
d.2 = d. = 0.202701 
d.3 = d3 = 0.0177127 



(A-6) 



(A-7) 



and 



3 

H’ = do + 2 X d = 2.8108034 (A-8) 

p=-3 

For a highpass filter with pass band from the cutoff frequency(B) to the 
maximum frequency, generate a low pass filter on the range 0 - B, and then 
subtract the central weight from unity and change the signs of the 
remainder of the weights. 

The filtered signal in the main program is generated by using the 
calcultaed filter weights as follows. 

do k=-m,m 

j=k 

if (k.lt.O) then 
j=k*(-l) 
endif 



(A-9) 

bb=-bk(j) 

if (k.eq.O) then 

bb=l.-bk(j) 

endif 

y(i)=y(i)+bb*s(i-k) 
end do 



5) hammg(dp,dt,pi,ain) 

This subroutine is a data tapering by using the modified Hamming 
window. It is often desirable to taper a random time series at each end to 
enhance certain characteristics of the spectral estimates. Tapering is 
multiplying the time series by a data window analogous to multiplying the 
correlation function by a lag window. Thus tapering the time series is 
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equivalent to applying a convolution operation to the raw Fourier 
transform. The purpose of tapering when viewed from its frequency 
domain effect is to suppress large side lobes in effective filter obtained with 
the raw transform. When looked at from the time domain, the object of 
tapering is to round off potential discontinuities at each end of the finite 
segment of the time history being analyzed. 

The used data window in this program is a modified Hamming 
window as given in Eq.(A-lO). 

0.54- 0.46 *cos(10jrt/D 0<t<T/10. 

W(t) = { 1.0 T/10 < t < 9T/10 (A-10) 

0.54 - 0.46 * cos( 10j:(T-t)/T) 9T/10 < t < T 



6) anal(dp,pi,ain,s) 

This subroutine converts a real signal to an analytic one by using the 
Hilbert transform given as follow, 



h(n) = - 



2 sin^(;m/2) 
K n 
0 , 



n 0, 
n = 0. 



(A-11) 



The Wigner distribution has a periodicity of N/2, where N is the 
number of data. Hence, even when the sampling of signal satisfies the 
Nyquist criteria, there are still aliasing components in the Wigner 
distribution function. If we sample the signal with Nyquist rate, its power 
density spectrum does not overlap with its own components. For Wigner 
distribution function, there are additional spectrum components causing 
the aliasing interference. However, a simple way to alleviate this problem 
for practical purpose is to increase the seperation of the spectrum groups by 
either doubling sampling rate or by interpolating with additional data 



points according to sampling theorem before the transform process. 
Another approach to avoid the aliasing is to use only the positive part the 
signal's frequency components, the analytical signal, before computing the 
Wigner distribution function. The analytic signal approach is convenient 
since it can be easily obtained with the Hilbert transform, which can also 
use the efficient FFT algorithm. 

This program uses the later case. The real time signal converts to 
complex variable by the Hilbert transform, that is, the imaginary part is 
generated by the convolution of the impulse response h(n) shown as follow, 

oo 

H{Sr(n)) = 5], h(n-m) Sr(m) (A-13) 

m=-oo 

where Sr is the original input signal. 

The analytic version of the real signal is made up of the real signal 
plus an imaginary part composed of the Hilbert transform of the real 
signal. 

s(t) = Sr(t) + j H{Sr(t)} (A-14) 



7) wigner(dp,dt,pi,s,c,mm,nn,wdf,mt,nf) 

This subroutine calculates the Wigner distribution function. The 
Wigner distribution function is obtained by using fast Fourier transform 
(FFT) of the local correlation of signal. At first, the local correlation of the 
signal is obtained and 2N points FFT is carried out. This program 
generates the results wdfTi,j) of the Wigner distribution function and 
pseudo Wigner-Ville distribution fimction on the given time series. The 
Wigner distribution has 1/2 of the frequency resolution of the power 
spectrum because 2N points FFT and the argument of the time signal and 
its conjugate contains a factor of 1/2. That is, df = l/(4*dp*dt). 
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To avoid the negative values and the elimination of the interference 
by the cross correlation term, this program uses the smoothing technique 
vidth Gaussian window function as following Eq.(A-15). 

G(t,w) = ! e (A-15) 

2 7T(7j (T^ 

The smoothing is obtained by the convolution integration the Wigner 
distribution and a Gaussian window function as given Eq(26). To perform 
the convolution on the sampled wdftij), the Gaussian function was first 
truncated so that it spanned the range ± 2(7j and ± 2(7^^. To avoid the 

negative value, Gaussian window size must have a larger or equal integer 
value mtand nf than (dp/7i)1^2 Por example, in the case of dp=1024, at least 
Gaussian window size must be 18x18 or larger. We recommend to select 
the square windows size. 



8) fft(dp,pi,c) 

This subroutine is a program about the Fast Fourier Transform by Jim 
Cooley's method. 
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APPENDIX B. PROGRAM LIST 



B-1. PROGRAM LIST ABOUT THE CALCULATION OF WDF 



PROGRAM PWVD 

***************************+********************************************* 



* 

* 

* 

* 



PSEUDO WIGNER-VILLE DISTRIBUTION FUNCTION 
PC VERSION 1.0 JAN. 1993 

G^y using highpass digital filter) 



* 

* 

* 



* * 
* 



* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

♦ 

* 

* 

* 

* 

* 

* 

* 



This is the program about the pseudo Wigner-Ville 
distribution(PWVD). The PWVD is a three dimensionaKtime, 
frequency, amplitude) representation of an input 
signal and is ideally suited for portraying transient 
phenomena. 



VARIABLES 



dp 



mm 

nn 

mt 

nf 



= the number of the sampled data point 
= output parameter for time 
= output parameter for fi*equency 
= smoothing parameter for time in Gaussian Fnc. 
= smoothing parameter for frequency 
in Gaussian Fnc. 

inname = input filename (sampled data) 
mvopt = option with respect to removing mean value 
if mvopt=l, then zero mean 
if mvopt=0, then no 
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* df = frequency resolution 

* dt = sampling interval 

* 

* ARRAYS 

* 



* 


tin(i) 


= sampled time data 


* 


ain(i) 


= sampled magnitude data 


* 


fain(i) 


= filtered data 


* 


wdf(i,j) 


= reduced array of PWVD 


* 


bk(i) 


= filter weights 


* 


s(i) 


= analytic signal 


* 


c(i) 


= local auto-correlation or the results of FFT 



* 

* * 

* VARIABLE DECLARATION * 

* * 

% ^ ^ ^ sjc ^ ^ nuttfiber of data points 

* Notes: For sample sizes greater than 2048, change 

* the parameter np. 

* 

integer dp,dp2,mvopt,redopt,mm,nn,nf,mt 
parameter (np=2048) 

real pi,tin(np),ain(np),wdf(256,128),dt,fain(np), 

1 bk(llOO) 

complex s(np*2),c(np*2) 
character*25 inname 
character as 
pi=atan(l.)*4. 

* — Print description of program — 

print* ’******** ****4=***4:*******:|:******* *+**♦♦*****• 

print*,' Pseudo Wigner-Ville Distribution’ 
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print* ’*******************************************' 

print* 

* — Set input parameters — 

print*,' Enter name of signal input file' 

read(5,901) inname 

print* 

print*, 'Number of the sampled data point' 
read(5,902) dp 

dp2=dp*2 

do 100 i=l,dp2 

100 s(i)=cmplx(0.,0.) 

♦ 

print*,' Do you wish to remove the mean value?' 
print*,' Enter 1 for yes or 0 for no' 
read(5,902) mvopt 
print* 

* 

print*. Do you want to apply a highpass digital' 
print*, 'filter to the original data ? (Y/N)' 
read(*,903) as 

if (as.eq.'Y'.or.as.eq.'y') then 

print*, 'Enter the cutoff frequency of 
print*, 'the digital highpass filter (in Hz) ' 
read(*,*) bw 

end 

if 

♦ 

fmin=0. 
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print*, Input the desired reduction size' 

print*,' input 1 for 64 by 32 ' 

print*,' input 2 for 128 by 64 ' 

print*,' input 3 for 128 by 128 ' 

print*,' input 4 for 256 by 128 ' 

read(5,902) redopt 

print* 

if (redopt.eq.l) then 
mm=dp/32 
nn=dp2/64 

elseif (redopt. eq.2) then 
mm = dp/64 
rm=dp2/128 

elseif (redopt. eq.3) then 
mm=dp/128 
nn=dp2/128 
else 

mm=dp/128 

nn=dp2/256 

endif 

print*, 'Input the desired smoothing window size 
print* 

print*,' input the smoothing parameter for frequency' 
print*,' in Gaussian function' 
read(5,902) nf 

print*,' input the smoothing parameter for time' 
print*,' in Gaussian function' 
read(5,902) mt 



* The calculation part of the program * 



Read the sampled data file 



♦ 



open(4, file=inname,status='old') 
call indata(dp,tin,ain) 
print*, 'finished subr. indata' 
close(4) 

* calculate the mean time interval 

call dtcalc(dp,tin,dt) 
print*, 'finished subr. dtcalc' 

if (mvopt.eq.l) then 
call mean(dp,ain) 
print*, 'finished subr. mean' 
endif 

♦ 

* Signal modifications 

* 

* Application of highpass digital filter 

if (as.eq.'Y'.or.as.eq.'y') then 
mo=dp/2 

* calculate the filter weighting 

call filter(mo,bw,dt,bk) 

* pass the highpass filter 

do 160 i=l,dp 
160 fain(i)=0. 

do 200 i=l,dp 
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do 170 k=-mo,mo 
j=k 

if (k.lt.O) then 
j=k*(-l) 
endif 

j=j+l 

ll=i-k 

if (ll.lt. 1) then 
ll=ll+dp 

elseif (ll.gt.dp) then 
ll=ll-dp 
endif 
bb=-bk(j) 
if (k.eq.O) then 
bb=l-bk(j) 
endif 

fain(i)=fain(i)+bb*ain(ll) 

170 continue 

200 continue 

do 210 i=l,dp 
210 ain(i)=fain(i) 

print*, 'finished digital filtering’ 
endif 

* Window application(modified hamming windows) 

call hammg(dp,dt,pi,ain) 
print*, 'finished subr. hammg' 

* Conversion of real signal to an analytic one 
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call anal(dp,pi,ain,s) 



print*, ’finished subr. anal' 



open(9,file='rwdf.out’,status='new’) 

open(10,file='rswdf.out',status=’new') 

* 

* Writing the WDF to output file * 

ttime=dp*dt 

df=l./(4.*dp*dt) 

fmax=2.*dp*df 

nx=dp2/nn 

ny=dp/mm 

write(9,904) tin(l),ttime,fmin,fmax 
write(9,*) nx,ny 

* Calculation of the Wigner distribution: 

print*, 'Calculating the PWVD' 
call wigner(dp,dt,pi,s,c,mm,nn,wdf,mt,nf) 
print*, 'finished wigner ville distribution function’ 
close(9) 

* 

* Writing of reduced & smoothed WVD to output file 

* 

write(10,904) tin(l),ttime,fmin,fmax 
writedO,*) nx,ny 
do 500 i=l,dp/mm 
do 500 j=l,dp2/nn 

write(10,905) wdRj,i) 

500 continue 

* Format statements 
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901 format(a25) 

902 format(i6) 

903 format(al) 

904 format(el2.5,/,el2.5/,el2.5/,el2.5) 

905 format(2x,el2.5) 

close(lO) 

print*, 'terminate the excution of Wigner-Ville distribution' 

stop 

end 



* * 

* SUBROUTINES * 

* * 



subroutine indata(dp,tin,ain) 

integer dp 
real tin(*),ain(*) 

********simple loop to read in time & amplitude************* 

do 100 j=l,dp 

read(4,*) tin(j) , ain(j) 

100 continue 
return 
end 






63 



subroutine dtcalc(dp,tin,dt) 

♦ 

* This subroutine calculates the delta t of the signal 

* 

* 

integer dp 
real tin(*),dt 

* 

dtsum = 0.0 

do 100 i = 1 , dp-1 

delt = tin(i-»-l) - tin(i) 
dtsum = dtsum -f- delt 

100 continue 

dt = dtsum / float(dp-l) 

return 

end 

:|c :<c :«c 4; 4: :)c :4c 4; 4; 4: 4; :)c ^ % :)c :<c :)c 4; :)c :)c :4c * :|c 4: :)c :«c :)c * 

subroutine mean(dp,ain) 

♦ 

* This subroutine calculates and removes the mean value 

* of the signal. 

* 



integer dp 
real ain(*),meanv 
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asum = 0.0 
do 100 i = 1 , dp 

asum = asum + ain(i) 
100 continue 

meanv = asum / dp 
do 200 i = 1 , dp 

ain(i) = ain(i) - meanv 
200 continue 

return 

end 



subroutine filter(mo,b,t,bk) 

* 

* Routine generates FIR filter weights. 

* Method devised by Potter, Bickford and Glaze. 

* There are a total of 2M+1 weights. ..filter generates M+1. 

* 

* — variables — 

* t = the sampling interval in second. 

* bw = cutofRhalf-power point) of the filter in Hz; 

* must be on the range from o to l/2t. 

* 

* Results are stored in bk 

* 

* — Note; in the case of highpass filter, the value of 

* weight bO must use 1-bO instead of bO. 

* 

dimension bk(*),d(3) 

data d0/0.355770 19/,d( l)/0.2436983/,d(2)/0.072 1 1497/, 

* d(3)/0.00630165/ 
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pi=atan(l.)*4. 

m=mo 

* first generate plain boxcar weights 
fact=2 *b*t 

bk(l)=fact 
fact=fact*pi 
do 5 i=l,m 
fi=i 

5 bk(i+l)=sin(fact*fi)/(pi*fi) 

* trapezoidal weighting at end 
bk(m+l)=bk(m+l)/2. 

* Now apply the Potter p3 10 window 
sumg=bk(l) 

do 15 i=l,m 
sum=d0 

fact=pi*float(i)/float(m) 
do 10 k=l,3 

10 sum=sum+2.*d(k)*cos(fact*float(k)) 

bk(i+ l)=bk(i+ l)*sum 
15 sumg=sumg+2.*bk(i+l) 

ml=m+l 
do 20 i=l,ml 
20 bk(i)=bk(i)/sumg 

return 
end 

subroutine hammg(dp,dt,pi,ain) 

* 

* This subroutine applies a modified hamming window 

* to the signal ain(t) 

* 

integer dp 
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real pi, ain(*),dt,mtime, del l,del2, const 

* 

mtime=(dp-l)*dt 

dell=0.1*mtime 

del2=0.9*mtime 

const=pi/dell 

do 100 j = 1 , dp 

t = (j-l)*dt 

if (t.le.dell) then 

ain(j) = ain(j) * (.54-0.46*cos(const*t)) 
elseif ((t.ge.del2).and.(t.le.mtime)) then 

ain(j)=ain(j)*(.54-0.46*cos(const*(mtime-t))) 

endif 

100 continue 

return 

end 

subroutine anal(dp,pi,ain,s) 

4c:4:5j::4::tc:tc5}c4c4c5ic5jc4t4c4c5fc:t:5f:5t:4c:4c5fc5}j5f::f:5lc5f:5ic5fc:t::4:5fc4:4:5lc:4;5icsjc:4c5lfijc:t;4:4:4:5fc5ic:^:5fc5jc:4c5fc4c5|c4c5i:s}c5lc:^c 

* 

* This subroutine converts a real signal to an 

* anal)d;ic one by using Hilbert transform. 

* 

* 



s(i) = analytic signal. 



integer dp 

real pi,ain(*),sum,sumb,val,sval 
complex s(*) 

do 100 i=l,dp 
sum=0.0 

do 200 j=l,dp 
sumb=0.0 



if (i-j.eq.O) go to 200 



n=i-j 

val=pi*n/2. 

sval=sin(val) 

sumb=ain(j)*sval*sval/val 

200 sum=sum+sumb 

s(i)=cmplx(ain(i),sum) 

100 continue 

return 

end 

subroutine wigner(dp,dt,pi,s,c,mm,nn,wdf,mt,nf) 

* 4::*: *:fc *****************:♦:**** :J:*:t:********5f:*:J:***** ♦♦****♦♦♦*♦ 

* 

* This subroutine calculates the WDF of the signal 

♦ 

* 
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integer dp,dp2 

real pi, dt,coef,wdf(256, 128), hg(-60:60, -60:60) 

complex s(*),dum,c(*) 

dp2 = dp* 2 

coef = 2.0 * dt 

df=l./(4.*dp*dt) 

nf2=nf*2 

mt2=mt*2 

fl=float(mt) 

f2=float(nf) 

* 

* Gaussian function 

* 

val=l./((2.*pi)**2*fl*f2*df*dt) 
do 20 j=-mt2,mt2 
ql=float(j) 
do 10 i=-nf2,nf2 
q2=float(i) 

Cf=-((ql*ql)/(2.*fl*fl))-((q2*q2)/(2.*f2*f2)) 
hg(i,j)=val*exp(cD 
10 continue 

20 continue 

* 

* initialize wdflij) 

* 

do 100 i=l,256 
do 100 j= 1,128 
100 wdflij)=0. 

do 5000 j = 1 , dp 

* local auto-correlation 

do 1000 i = 1 , dp+1 
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1000 



1400 

1500 



if 0-ge.i) then 
dum = s(j-i+l) 
else 

dum = cmplx(0.,0.) 
endif 

c(i) = coef * (s(j+i-l)*conjg(dum)) 
if (i.ne.l.and.i.ne.dp+1) then 
c(dp2-i+2) = conjg(c(i)) 
endif 
continue 

call fft(dp,pi,c) 

ik=mod((j-l),mm) 
if (ik.eq.O) then 

do 1500 i=l,dp2,nn 

write(9,1400) real(c(i)) 
format(2x,el2.5) 
continue 
endif 
ml=0 
do 4000 
m=l,dp,mm 
ml=ml+l 
nl=0 

if (abs(j-m).le.mt2) then 
do 3000 n=l,dp2,nn 
nl=nl+l 

do 2500 kk=n-nf2,n+nf2 
kl=kk 

iftkk.lt. 1) kl=kk+dp2 
iftkk.gt.dp2) kl=kk-dp2 

wdftnl,ml)=wdftnl,ml)+real(c(kl))*hg(kk-nj-m)*df*dt 
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2500 


continue 


3000 


continue 




endif 


4000 


continue 


5000 


continue 



return 

end 



subroutine fTt(dp,pi,c) 

* 

* This subroutine is the Fast Fourier Transform 

* 

integer dp,dp2,val,coef,coefl 
real pi 

complex dum,c(*),dum3,dum2 

* 

dp2 = dp*2 
const=float(dp2) 
val=alog(const)/alog(2.)+. 1 

j=l 

do 40 i=l,dp2-l 
if (i.ge.j) go to 10 
dum3=c(j) 
c(j)=c(i) 
c(i)=dum3 
10 k=dp 

20 if (k.ge.j) go to 30 
j=j-k 
k=k/2 
go to 20 
30 j=j+k 
40 continue 



do 70 n=l,val 
coef=2**n 
coefl=coef/2 
dum2=cmplx(l.,0.) 
theta=pi/float(coefl) 
dum=cmplx(cos(theta),-sin(theta)) 
do 60 j=l,coefl 
do 50 i=j,dp2,coef 
ii=i+coefl 
dum3=c(ii)*du 
m2 c(ii)=c(i)- 
dum3 

c(i)=c(i)+dum3 



50 


continue 




dum2=dum2 


60 


continue 


70 


continue 




return 




end 



72 



B-2. PROGRAM LIST OF 3-D PLOT 



PROGRAM 3DPLOT 

This program uses the graphic package CA-DISSPLA to plot 
the results of WDF 



* 


ttime 


= time record length 


* 


tini 


= initial time 


* 


fmin 


= start frequency 


* 


fmax 


= stop frequency 


* 


nx 


= the number of the frequency data at time n x dt 


* 


ny 


= the number of the time data at frequency m x df 



* 



* Declaring variables 

* 

real rwdf(32768) 
integer nx,ny 
character*25 fname 

write(*,*) 'input file name ?' 
read(*,20) fname 

20 format(a25) 

* 

open(15,file=fname,status='old') 

read(15,*) tini 
read(15,*) ttime 
read(15,*) fmin 
read(15,*) fmax 
read(15,*) nx,ny 

write(*,*) tini, ttime, fmin, fmax, nx,ny 
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n=nx*ny 
do 100 i=l,n 

read(15,*) rwdfli) 

100 continue 

close(15) 

smax=rwdf(l) 
do 200 i=l,n 

if (smax.lt.rwdfli)) then 
smax=rwdf(i) 
endif 

200 continue 

write(*,*) 'smax =’, smax 

writeC*,*) 'input maximum z-axis value’ 

read(*,*) fac 

call pdev('ln03', ieer) 

* plotting 

call hwshd 
call swissm 

call shdchKOO., 1,0.002,1) 
call height(0.2) 
call physor{0.7,0.625) 
call area2d(7.5,9.75) 

call messagCWIGNER-VILLE DISTRIBUTION $’,100,1.1,8.2) 
call blsur 

call volm3d(8.,8.,9.) 

call x3name(’Frequency (Hz) $’,100) 

call y3name(’Time (sec) $’,100) 
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call z3name(’Amplitude $’,100) 

call vuangl(-60.,30.,30.) 

call zaxangOO.) 

tstep=ttime/4. 

fstep=(fmax-fmin)/4. 

tmax=tini+ttime 

call graf3d(fimn,fstep,fmax,tini,tstep,tmax,0.,'SCALE',fac) 

call surmat(rwdf,l,nx,l,ny,l) 

call end3gr(0) 
call endpKo) 
call donepl 

stop 

end 
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B-3. PROGRAM LIST OF CONTOUR PLOT 



PROGRAM CONTOUR 

This program uses the graphic package CA-DISSPLA to plot 
the results of WDF 



♦ 


ttime 


= time record length 


* 


tini 


= initial time 


* 


fmin 


= start frequency 


* 


fmax 


= stop frequency 


* 


nx 


= the number of the frequency data at time n x dt 


* 


ny 


= the number of the time data at fi*equency m x df 



* 



* Declaring variables 

* 

real r(64,128) 
common work(9000) 
integer nx,ny 
character*25 fname 

write(*,*) 'input file name ?' 
read(*,20) fname 

20 format(a25) 

* 

open(15,file=fname,status=’old') 

read(15,*) tini 
read(15,*) ttime 
read(15,*) fmin 
read(15,*) fmax 
read(15,*) nx,ny 

write(*,*) tini, ttime, fmin, fmax, nx,ny 
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smax=0. 
do 100 i=l,ny 

do 100 j=l,nx 
read(15,*) r(ij) 
if (r(iJ).lt.O.) r(ij)=0. 
if (smax.lt. r(ij)) then 
smax=r(i,j) 
endif 

100 continue 

close(15) 

* Normalizing 

do 200 i=l,ny 
do 200 j=l,nx 

r(iJ)=r(ij)/smax 
200 continue 

write(*,*) 'smax =', smax 

write(*,*) 'input maximum z-axis value' 

read(*,*) fac 

call pdev('ln03', ieer) 

* plotting 

call hwshd 
call swissm 

call shdchr(90., 1,0.002,1) 
call height(0.2) 
call page(8.5,ll.) 
call physoKl.5,1.5) 
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call area2d(6.,6.) 

call headinCWIGNER-VILLE DISTRIBUTION CONTOURS', 
* 100 , 1 . 1 , 1 ) 
call blsur 

call ynameCFrequency (Hz) $’,100) 

call xnameCTime (sec) $’,100) 

call yaxang(O-) 

tstep=ttime/5. 

fstep=(fmax-fmin)/5. 

tmax=tini+ttime 

call graf(tini,tstep,tmax,fimn,fstep,fmax) 

call frame 

call bcomonOOOO) 

scale=l./30. 

call conmak(r,ny,nx, scale) 
contour plot 

call conlin(0,’SOLID’, ’NOLABELS’,1,3) 
call conlin(l.,’DASH’, ’NOLABELS’, 1,1) 
call conang(50.) 
call raspln(0.25) 

call contur(2,’LABELS’,’DRAW’) 

call endpKo) 
call donepl 

stop 

end 
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B-4. PROGRAM LIST OF PLOT OF TIME SERIES 



PROGRAM TIMEPLOT 

* 

* This program uses the graphic package CA-DISSPLA to plot 

* the time series. 

* 

* Declaring variables 

* 

dimension x(4000),y(4000) 
character*60 title 
character*25 inname 

write(*,*) 'input file name ?' 
read(*,1000) inname 

1000 format(a25) 
write(*,*) 'title ?' 
read(*,1001) title 

1001 format(a60) 

write(*,*) 'the number of the data point' 
read(*,*) np 

write(*,*) 'maximum scale of the time in a figure' 
read(*,*) xmax 

write(*,*) 'maximum scale of the magnitude in the figure' 
read(*,*) ymax 

* 

open(8,file=inname,status='old') 

do 100 i=l,np 

read(8,*) x(i),y(i) 

100 continue 
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close(8) 



call pdev(’ln03’, ieer) 

plotting 

call hwshd 
call swissm 

call shdchr(90., 1,0.002,1) 

call height(0.2) 

call page(8.5,ll.) 

call physor(1.5,1.5) 

call area2d(6.,6.) 

call xnameCTime (sec) $’,100) 

call ynameCFrequency (Hz) $’,100) 

call headin(title,60,l.l,l) 

call thkfrm(O.Ol) 

call yaxangOO.) 

call graf(0.,xmax/4.,xmax,-ymax,0.5,ymax) 

call grid(l,l) 

call curve(x,y,np,0) 

call endpl(o) 
call donepl 

stop 

end 
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